In [1]:
import astropy.io.fits as fits
import numpy as np
from glob import glob
from matplotlib import pyplot as plt
import matplotlib.colors as colors
from astropy.visualization import simple_norm
import ccdproc as ccd
from astropy import visualization as viz
from astropy.visualization import (imshow_norm, MinMaxInterval, SqrtStretch, ZScaleInterval)
from astropy.coordinates import SkyCoord
from photutils.aperture import ApertureStats
from astropy.stats import SigmaClip
from photutils.aperture import aperture_photometry, CircularAperture, CircularAnnulus
from photutils.detection import DAOStarFinder
from astropy.stats import mad_std
from astropy.wcs import WCS
from astropy.table import Table
from astropy.io import ascii
from astropy import units as u
from datetime import datetime
from astropy.time import Time
from astropy.modeling import models, fitting
from scipy.optimize import curve_fit
from astropy.visualization import ImageNormalize
from astropy.coordinates import solar_system_ephemeris, EarthLocation
from astropy.coordinates import get_body_barycentric, get_body, get_moon

Data Reduction¶

Files¶

In [2]:
directory = '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/*'
glob(directory)
Out[2]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/alberio',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/bias',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/calibration',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/dark',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/deneb',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/epsPeg',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/europa',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/flat',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/jupiter',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/sun',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/vega']
In [3]:
directory2 = '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/mygroup/*'
glob(directory2)
Out[3]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/mygroup/bias',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/mygroup/calibration',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/mygroup/dark',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/mygroup/flat',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/mygroup/sky',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/mygroup/vega']
In [4]:
files_bias = glob(glob(directory)[1]+'/*')
files_dark = glob(glob(directory)[3]+'/*')
files_flat = glob(glob(directory)[7]+'/*')
files_cali = glob(glob(directory)[2]+'/*')
files_albr = glob(glob(directory)[0]+'/*')
files_denb = glob(glob(directory)[4]+'/*')
files_epsP = glob(glob(directory)[5]+'/*')
files_euro = glob(glob(directory)[6]+'/*')
files_jupi = glob(glob(directory)[8]+'/*')
files_sunn = glob(glob(directory)[9]+'/*')
files_vega = glob(glob(directory)[10]+'/*')
In [5]:
files_bias2 = glob(glob(directory2)[0]+'/*')
files_dark2 = glob(glob(directory2)[2]+'/*')
files_flat2 = glob(glob(directory2)[3]+'/*')
files_cali2 = glob(glob(directory2)[1]+'/*')
In [6]:
ccd_shape = fits.getdata(files_bias[0]).shape
In [ ]:
In [ ]:
In [ ]:

Bias¶

In [7]:
bias = np.zeros((len(files_bias), ccd_shape[0], ccd_shape[1]))
for file in range(len(files_bias)):
    bias[file] = fits.getdata(files_bias[file])
In [8]:
master_bias = np.median(bias, axis=0)
imshow_norm(master_bias, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[8]:
(<matplotlib.image.AxesImage at 0x7f2ab4d3a580>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2abd18ad00>)
In [9]:
bias2 = np.zeros((len(files_bias2), ccd_shape[0], ccd_shape[1]))
for file in range(len(files_bias2)):
    bias2[file] = fits.getdata(files_bias2[file])
In [10]:
master_bias2 = np.median(bias2, axis=0)
imshow_norm(master_bias2, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[10]:
(<matplotlib.image.AxesImage at 0x7f2ab343d7c0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab4c8dbb0>)

Dark¶

Other Group¶

In [11]:
files_dark_07 = glob(files_dark[2]+'/*')
files_dark_15 = glob(files_dark[0]+'/*')
files_dark_30 = glob(files_dark[3]+'/*')
files_dark_40 = glob(files_dark[1]+'/*')
In [12]:
dark_07 = np.zeros((len(files_dark_07), ccd_shape[0], ccd_shape[1]))
exptime_dark_07 = fits.getheader(files_dark_07[0])['EXPTIME']
for file in range(len(files_dark_07)):
    dark_07[file] = fits.getdata(files_dark_07[file])
In [13]:
master_dark_07 = np.median(dark_07-master_bias, axis=0)
imshow_norm(master_dark_07, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[13]:
(<matplotlib.image.AxesImage at 0x7f2ab33bd2b0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab33fa880>)
In [14]:
dark_15 = np.zeros((len(files_dark_15), ccd_shape[0], ccd_shape[1]))
exptime_dark_15 = fits.getheader(files_dark_15[0])['EXPTIME']
for file in range(len(files_dark_15)):
    dark_15[file] = fits.getdata(files_dark_15[file])
In [15]:
master_dark_15 = np.median(dark_15-master_bias, axis=0)
imshow_norm(master_dark_15, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[15]:
(<matplotlib.image.AxesImage at 0x7f2ab33a73a0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab34622e0>)
In [16]:
dark_30 = np.zeros((len(files_dark_30), ccd_shape[0], ccd_shape[1]))
exptime_dark_30 = fits.getheader(files_dark_30[0])['EXPTIME']
for file in range(len(files_dark_30)):
    dark_30[file] = fits.getdata(files_dark_30[file])
In [17]:
master_dark_30 = np.median(dark_30-master_bias, axis=0)
imshow_norm(master_dark_30, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[17]:
(<matplotlib.image.AxesImage at 0x7f2ab330e580>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab333a910>)
In [18]:
dark_40 = np.zeros((len(files_dark_40), ccd_shape[0], ccd_shape[1]))
exptime_dark_40 = fits.getheader(files_dark_40[0])['EXPTIME']
for file in range(len(files_dark_40)):
    dark_40[file] = fits.getdata(files_dark_40[file])
In [19]:
master_dark_40 = np.median(dark_40-master_bias, axis=0)
imshow_norm(master_dark_40, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[19]:
(<matplotlib.image.AxesImage at 0x7f2ab327a760>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab3320f70>)

My Group¶

In [20]:
files_dark2_10 = glob(files_dark2[3]+'/*')
files_dark2_20 = glob(files_dark2[4]+'/*')
files_dark2_40 = glob(files_dark2[1]+'/*')
files_dark2_50 = glob(files_dark2[0]+'/*')
files_dark2_60 = glob(files_dark2[2]+'/*')
In [21]:
dark2_10 = np.zeros((len(files_dark2_10), ccd_shape[0], ccd_shape[1]))
exptime_dark2_10 = fits.getheader(files_dark2_10[0])['EXPTIME']
for file in range(len(files_dark2_10)):
    dark2_10[file] = fits.getdata(files_dark2_10[file])
In [22]:
master_dark2_10 = np.median(dark2_10-master_bias2, axis=0)
imshow_norm(master_dark2_10, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[22]:
(<matplotlib.image.AxesImage at 0x7f2ab32705e0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab32c19d0>)
In [23]:
dark2_20 = np.zeros((len(files_dark2_20), ccd_shape[0], ccd_shape[1]))
exptime_dark2_20 = fits.getheader(files_dark2_20[0])['EXPTIME']
for file in range(len(files_dark2_20)):
    dark2_20[file] = fits.getdata(files_dark2_20[file])
In [24]:
master_dark2_20 = np.median(dark2_20-master_bias2, axis=0)
imshow_norm(master_dark2_20, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[24]:
(<matplotlib.image.AxesImage at 0x7f2ab023c7f0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab027ddf0>)
In [25]:
dark2_40 = np.zeros((len(files_dark2_40), ccd_shape[0], ccd_shape[1]))
exptime_dark2_40 = fits.getheader(files_dark2_40[0])['EXPTIME']
for file in range(len(files_dark2_40)):
    dark2_40[file] = fits.getdata(files_dark2_40[file])
In [26]:
master_dark2_40 = np.median(dark2_40-master_bias2, axis=0)
imshow_norm(master_dark2_40, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[26]:
(<matplotlib.image.AxesImage at 0x7f2ab01a17f0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab3340f40>)
In [27]:
dark2_50 = np.zeros((len(files_dark2_50), ccd_shape[0], ccd_shape[1]))
exptime_dark2_50 = fits.getheader(files_dark2_50[0])['EXPTIME']
for file in range(len(files_dark2_50)):
    dark2_50[file] = fits.getdata(files_dark2_50[file])
In [28]:
master_dark2_50 = np.median(dark2_50-master_bias2, axis=0)
imshow_norm(master_dark2_50, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[28]:
(<matplotlib.image.AxesImage at 0x7f2ab018a9a0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab024a2e0>)
In [29]:
dark2_60 = np.zeros((len(files_dark2_60), ccd_shape[0], ccd_shape[1]))
exptime_dark2_60 = fits.getheader(files_dark2_60[0])['EXPTIME']
for file in range(len(files_dark2_60)):
    dark2_60[file] = fits.getdata(files_dark2_60[file])
In [30]:
master_dark2_60 = np.median(dark2_60-master_bias2, axis=0)
imshow_norm(master_dark2_60, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[30]:
(<matplotlib.image.AxesImage at 0x7f2ab00f4b50>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab01ecfa0>)

Flat¶

Other Group¶

In [ ]:
In [31]:
flat = np.zeros((len(files_flat), ccd_shape[0], ccd_shape[1]))
exptime_flat = fits.getheader(files_flat[0])['EXPTIME']
for file in range(len(files_flat)):
    flat[file] = fits.getdata(files_flat[file])
In [32]:
master_flat = np.median(flat-master_bias-(exptime_flat/exptime_dark_07)*master_dark_07, axis=0)
imshow_norm(master_flat, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[32]:
(<matplotlib.image.AxesImage at 0x7f2ab0060e50>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab00f4ee0>)
In [33]:
np.mean(master_flat)
Out[33]:
22.085019131835743
In [34]:
master_norm_flat = master_flat / np.mean(master_flat)
np.mean(master_norm_flat)
Out[34]:
0.9999999999999992

My Group¶

In [35]:
flat2 = np.zeros((len(files_flat2), ccd_shape[0], ccd_shape[1]))
exptime_flat2 = fits.getheader(files_flat2[0])['EXPTIME']
for file in range(len(files_flat2)):
    flat2[file] = fits.getdata(files_flat2[file])
In [36]:
master_flat2 = np.median(flat2-master_bias2-(exptime_flat2/exptime_dark2_10)*master_dark2_10, axis=0)
imshow_norm(master_flat2, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[36]:
(<matplotlib.image.AxesImage at 0x7f2aaffd5580>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ab0018280>)
In [37]:
master_norm_flat2 = master_flat2 / np.mean(master_flat2)
np.mean(master_norm_flat2)
Out[37]:
1.0000000000000004

Wavelength Calibration¶

$[master\_flat-\frac{exptime\_flat}{exptime\_dark}(master\_dark)]_{ave}[\frac{(object-master\_bias) -\frac{exptime\_obj}{exptime\_dark}(master\_dark)}{(master\_flat)-\frac{exptime\_flat}{exptime\_dark}(master\_dark)}]_{pix}$

$[(flat-bias)-\frac{t_f}{t_d}(dark-bias)]_{ave}[\frac{(object-bias)-\frac{t_o}{t_d}(dark-bias)}{(flat-bias)-\frac{t_f}{t_d}(dark-bias)}]_{pix}$

In [38]:
def lin_fit(x, m, b):
    return m*x+b
def quad_fit(x, a, b, c):
    return a*x**2 + b*x + c
def gauss_fit(x, a, x0, sigma):
    return a*np.exp(-1*(x-x0)**2/(2*sigma**2))

Other Group¶

In [39]:
files_cali
Out[39]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/calibration/cal_Hg_7s.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/calibration/cal_neon_10s.FIT']
In [40]:
cali_Hg = fits.getdata(files_cali[0])
In [41]:
exptime_Hg = fits.getheader(files_cali[0])['EXPTIME']

final_Hg = (cali_Hg - master_bias - (exptime_Hg/exptime_dark_07)*master_dark_07)/(master_norm_flat)

/tmp/ipykernel_5556/1387433607.py:3: RuntimeWarning: divide by zero encountered in true_divide
  final_Hg = (cali_Hg - master_bias - (exptime_Hg/exptime_dark_07)*master_dark_07)/(master_norm_flat)
/tmp/ipykernel_5556/1387433607.py:3: RuntimeWarning: invalid value encountered in true_divide
  final_Hg = (cali_Hg - master_bias - (exptime_Hg/exptime_dark_07)*master_dark_07)/(master_norm_flat)
In [42]:
imshow_norm(final_Hg[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[42]:
(<matplotlib.image.AxesImage at 0x7f2aaffbbfd0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2aaff95220>)
In [43]:
em_Hg = np.nanmedian(final_Hg[375:425,:], axis=0)
px_Hg = np.arange(len(em_Hg))
em_Hg[em_Hg==np.nan] = 0
plt.plot(em_Hg)
# 435.835 nm (blue), 546.074 nm (green)
Out[43]:
[<matplotlib.lines.Line2D at 0x7f2aaff1fca0>]
In [44]:
#plt.plot(em_Hg[0:26])
#plt.plot(em_Hg[26:50])
#plt.plot(em_Hg[160:200])
plt.plot(em_Hg[300:360])
#plt.plot(em_Hg[500:600])
#plt.plot(em_Hg[690:730])
#plt.plot(em_Hg[825:900])
Out[44]:
[<matplotlib.lines.Line2D at 0x7f2aafe9e550>]
In [45]:
cali_Ne = fits.getdata(files_cali[1])
In [46]:
exptime_Ne = fits.getheader(files_cali[1])['EXPTIME']

final_Ne = (cali_Ne - master_bias - (exptime_Ne/exptime_dark_07)*master_dark_07)/(master_norm_flat)

/tmp/ipykernel_5556/1325804219.py:3: RuntimeWarning: divide by zero encountered in true_divide
  final_Ne = (cali_Ne - master_bias - (exptime_Ne/exptime_dark_07)*master_dark_07)/(master_norm_flat)
/tmp/ipykernel_5556/1325804219.py:3: RuntimeWarning: invalid value encountered in true_divide
  final_Ne = (cali_Ne - master_bias - (exptime_Ne/exptime_dark_07)*master_dark_07)/(master_norm_flat)
In [47]:
imshow_norm(final_Ne[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[47]:
(<matplotlib.image.AxesImage at 0x7f2aafe7bee0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2aafed0790>)
In [48]:
em_Ne = np.nanmedian(final_Ne[375:425,:], axis=0)
plt.plot(em_Ne)
Out[48]:
[<matplotlib.lines.Line2D at 0x7f2aafde2760>]
In [49]:
plt.plot(em_Ne[:400])
#plt.plot(em_Ne[0:20])
#plt.plot(em_Ne[20:50])
#plt.plot(em_Ne[190:230])
#plt.plot(em_Ne[230:250])
#plt.plot(em_Ne[330:355])
#plt.plot(em_Ne[355:400])
Out[49]:
[<matplotlib.lines.Line2D at 0x7f2aafdc4ca0>]

My Group¶

In [ ]:
In [50]:
files_cali2
Out[50]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/mygroup/calibration/Hydrogen_10s.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/mygroup/calibration/Mercury_20s.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/mygroup/calibration/Neon_10s.FIT']
In [51]:
cali2_H  = fits.getdata(files_cali2[0])
cali2_Hg = fits.getdata(files_cali2[1])
cali2_Ne = fits.getdata(files_cali2[2])
In [52]:
gain = fits.getheader(files_cali2[1])['EGAIN']
In [53]:
exptime2_Hg = fits.getheader(files_cali2[1])['EXPTIME']
final2_Hg = (cali2_Hg - master_bias2 - (exptime2_Hg/exptime_dark2_20)*master_dark2_20)/(master_norm_flat2)

/tmp/ipykernel_5556/786135649.py:2: RuntimeWarning: divide by zero encountered in true_divide
  final2_Hg = (cali2_Hg - master_bias2 - (exptime2_Hg/exptime_dark2_20)*master_dark2_20)/(master_norm_flat2)
/tmp/ipykernel_5556/786135649.py:2: RuntimeWarning: invalid value encountered in true_divide
  final2_Hg = (cali2_Hg - master_bias2 - (exptime2_Hg/exptime_dark2_20)*master_dark2_20)/(master_norm_flat2)
In [54]:
imshow_norm(final2_Hg[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[54]:
(<matplotlib.image.AxesImage at 0x7f2aafd3b160>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2aafd81610>)
In [55]:
em2_Hg = np.nanmedian(final2_Hg[375:425,:], axis=0)*gain/exptime2_Hg
px2_Hg = np.arange(len(em2_Hg))
plt.plot(em2_Hg)
Out[55]:
[<matplotlib.lines.Line2D at 0x7f2aafd0ecd0>]
In [56]:
exptime2_H = fits.getheader(files_cali2[0])['EXPTIME']
final2_H = (cali2_H - master_bias2 - (exptime2_H/exptime_dark2_10)*master_dark2_10)/(master_norm_flat2)

/tmp/ipykernel_5556/2111062723.py:2: RuntimeWarning: divide by zero encountered in true_divide
  final2_H = (cali2_H - master_bias2 - (exptime2_H/exptime_dark2_10)*master_dark2_10)/(master_norm_flat2)
/tmp/ipykernel_5556/2111062723.py:2: RuntimeWarning: invalid value encountered in true_divide
  final2_H = (cali2_H - master_bias2 - (exptime2_H/exptime_dark2_10)*master_dark2_10)/(master_norm_flat2)
In [57]:
imshow_norm(final2_H[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[57]:
(<matplotlib.image.AxesImage at 0x7f2aafc81cd0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2aafcca7c0>)
In [58]:
em2_H = np.nanmedian(final2_H[375:425,:], axis=0)*gain/exptime2_H
px2_H = np.arange(len(em2_H))
plt.plot(em2_H)
Out[58]:
[<matplotlib.lines.Line2D at 0x7f2aafbe4880>]
In [ ]:
In [59]:
exptime2_Ne = fits.getheader(files_cali2[1])['EXPTIME']
final2_Ne = (cali2_Ne - master_bias2 - (exptime2_Ne/exptime_dark2_10)*master_dark2_10)/(master_norm_flat2)

/tmp/ipykernel_5556/300449279.py:2: RuntimeWarning: divide by zero encountered in true_divide
  final2_Ne = (cali2_Ne - master_bias2 - (exptime2_Ne/exptime_dark2_10)*master_dark2_10)/(master_norm_flat2)
/tmp/ipykernel_5556/300449279.py:2: RuntimeWarning: invalid value encountered in true_divide
  final2_Ne = (cali2_Ne - master_bias2 - (exptime2_Ne/exptime_dark2_10)*master_dark2_10)/(master_norm_flat2)
In [60]:
imshow_norm(final2_Ne[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[60]:
(<matplotlib.image.AxesImage at 0x7f2aafbd1520>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2aafc09f10>)
In [61]:
em2_Ne = np.nanmedian(final2_Ne[375:425,:], axis=0)*gain/exptime2_Ne
px2_Ne = np.arange(len(em2_Ne))
plt.plot(em2_Ne)
Out[61]:
[<matplotlib.lines.Line2D at 0x7f2aafb37070>]
In [62]:
#plt.plot(em2_Hg, label='Hg')
#plt.plot(em2_Ne, label='Ne')
#plt.plot(wavelengths, em2_Hg, label='Hg')
#plt.plot(em2_Ne, label='Ne')
#plt.plot(em2_H, label='H')
#plt.ylim(top=50, bottom=-10)
#plt.xlim(left=395, right=415)
#plt.legend()

#There are two other blue lines at 404.656 nm and 407.781 nm and a weak line at 491.604 nm.
# The prominent mercury lines are at 435.835 nm (blue), 546.074 nm (green), and a pair at 576.959 nm and 579.065 nm (yellow-orange).
In [63]:
px_H_633 = 434.047#*u.nm # blue
px_H_382 = 486.133#*u.nm # teal
px_Hg_622 = 435.835#*u.nm # blue
px_Hg_93  = 546.074 #*u.nm # green
px_Ne_121 = 540.05616# *u.nm
px_Ne_151 = 534.10938 #*u.nm

#plt.plot(em2_H[:])
#plt.plot(em2_Hg[:])
#plt.plot(em2_Ne[:])
#plt.ylim(top=200, bottom=-5)

#np.argmax(em2_Ne[140:190])

# The prominent mercury lines are at 435.835 nm (blue), 546.074 nm (green), 
# and a pair at 576.959 nm and 579.065 nm (yellow-orange). There are two other 
# blue lines at 404.656 nm and 407.781 nm and a weak line at 491.604 nm
In [64]:
x = [633, 382, 622, 93, 121, 151]
wave = [px_H_633, px_H_382, px_Hg_622, px_Hg_93, px_Ne_121, px_Ne_151]
fit = curve_fit(lin_fit, x, wave)[0]
In [65]:
print(fit[0], fit[1])

-0.2078095090608855 565.3815295233154
In [66]:
ele = ['Hydrogen', 'Hydrogen', 'Mercury', 'Mercury', 'Neon', 'Neon']
names = ['Pixel', 'Wavelength', 'Element']
tblr_fit = Table([x*u.pixel, wave*u.nm, ele], names=names)
tblr_fit.write('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/table_wavelengthsolu.txt', format='latex', overwrite=True)
In [67]:
pwd
Out[67]:
'/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix'
In [68]:
px = np.arange(len(em2_Hg))
wavelengths = lin_fit(px, fit[0], fit[1])
In [69]:
plt.plot(px, wavelengths, '--', color='red')
plt.scatter(x, wave, color='black')
plt.xlim(left=0, right=800)
plt.ylim(bottom=350)
plt.xlabel('Pixel Position (pix)')
plt.ylabel('Wavelength (nm)')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/wavevpx.pdf', bbox_inches='tight')
In [70]:
fit_c = curve_fit(quad_fit, x, wave)[0]
wavelengths_c = quad_fit(px, fit_c[0], fit_c[1], fit_c[2])
In [71]:
fit_c[0]
Out[71]:
-2.2193685731419777e-06
In [72]:
np.max(wavelengths)
Out[72]:
565.3815295233154
In [73]:
plt.plot(wavelengths, em2_H)
Out[73]:
[<matplotlib.lines.Line2D at 0x7f2ac2315b20>]
In [74]:
plt.plot(wavelengths_c, em2_H)
Out[74]:
[<matplotlib.lines.Line2D at 0x7f2ac22cdbb0>]
In [75]:
px_H_633 = 434.047#*u.nm # blue
px_H_382 = 486.133#*u.nm # teal
px_Hg_622 = 435.835#*u.nm # blue
px_Hg_93  = 546.074 #*u.nm # green
px_Ne_121 = 540.05616# *u.nm
px_Ne_151 = 534.10938 #*u.nm
H_balmer62 = 410.1734
H_balmer72 = 397.0075
He_line = 537.02992
H_balmer82 = 388.9064

fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em2_H, label='Hydrogen', color='tab:blue')
plt.plot(wavelengths, em2_Hg, label='Mercury', color='tab:orange')
plt.plot(wavelengths, em2_Ne, label='Neon', color='tab:green')
plt.title('Spectrum of Calibrators')
plt.xlabel('Wavelength (nm)')
plt.ylabel('Intensity')
plt.text(px_H_633-7, 50, 'Balmer (5-2)', rotation='vertical', color='tab:blue')
plt.text(px_H_382+2, 150, 'Balmer (4-2)', rotation='vertical', color='tab:blue')
plt.text(px_Hg_622+2, 150, 'Hg I (62-39)', rotation='vertical', color='tab:orange')
plt.text(px_Hg_93+2, 50, 'Hg I (62-44)', rotation='vertical', color='tab:orange')
plt.text(px_Ne_121, 20, 'Ne I (152-134)', rotation='vertical', color='tab:green')
plt.text(px_Ne_151, 10, 'Ne I (166-148)', rotation='vertical', color='tab:green')
plt.legend()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/calibrators_spectrum.pdf', bbox_inches='tight')
In [76]:
fig = plt.figure()
gs = fig.add_gridspec(3, hspace=0)
axs = gs.subplots(sharex=True, sharey=False)
fig.set_figheight(4)
fig.set_figwidth(10)
axs[0].set_xlim(left=np.min(wavelengths), right=np.max(wavelengths))

axs[0].plot(wavelengths, em2_H, label='Hydrogen', color='black')
axs[1].plot(wavelengths, em2_Hg, label='Mercury', color='black')
axs[2].plot(wavelengths, em2_Ne, label='Neon', color='black')

axs[0].text(px_H_633+2,  50, 'H γ', rotation='vertical', color='tab:blue')
axs[0].text(px_H_382+2,  50, 'H β', rotation='vertical', color='tab:blue')
axs[1].text(px_Hg_622+2, 50, 'Hg I', rotation='vertical', color='tab:orange')
axs[1].text(px_Hg_93+2,  50, 'Hg I', rotation='vertical', color='tab:orange')
axs[2].text(px_Ne_121+2, 5, 'Ne I', rotation='vertical', color='tab:green')
axs[2].text(px_Ne_151, 5, 'Ne I', rotation='vertical', color='tab:green')

for ax in fig.get_axes():
    ax.label_outer()
    ax.legend(loc='upper left')
    
fig.supylabel('Flux (e-/s)')
fig.supxlabel('Wavelength (nm)')    

    
fig.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/calibrators_spectrum.pdf', bbox_inches='tight')
In [77]:
fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em2_H, label='Hydrogen', color='black')
plt.text(px_H_633+3, 50, 'Balmer (5-2)', rotation='vertical')
plt.text(px_H_382+3, 150, 'Balmer (4-2)', rotation='vertical')
plt.title('Spectrum of Hydrogen')
plt.ylabel('Flux (e-/s)')
plt.xlabel('Wavelength (nm)')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/hydrogen_spectrum.pdf', bbox_inches='tight')
In [78]:
fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em2_Hg, label='Mercury', color='black')
plt.text(px_Hg_622+3, 100, 'Hg I (62-39)', rotation='vertical')
plt.text(px_Hg_93+3, 50, 'Hg I (62-44)', rotation='vertical')
plt.title('Spectrum of Mercury')
plt.ylabel('Flux (e-/s)')
plt.xlabel('Wavelength (nm)')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/mercury_spectrum.pdf', bbox_inches='tight')
In [79]:
fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em2_Ne, label='Neon', color='black')
plt.text(px_Ne_121+3, 6, 'Ne I (152-134)', rotation='vertical')
plt.text(px_Ne_151-3, 4, 'Ne I (166-148)', rotation='vertical')
plt.title('Spectrum of Neon')
plt.ylabel('Flux (e-/s)')
plt.xlabel('Wavelength (nm)')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/neon_spectrum.pdf', bbox_inches='tight')

Science¶

Data¶

In [80]:
exptime_alberio = fits.getheader(files_albr[0])['EXPTIME']
files_albr
Out[80]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/alberio/alberioa_40s_1.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/alberio/alberioa_40s_2.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/alberio/alberioa_40s_3.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/alberio/alberioa_40s_4.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/alberio/alberioa_40s_5.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/alberio/alberiob_40s_1.FIT']
In [81]:
data_alb = np.zeros((len(files_albr), ccd_shape[0], ccd_shape[1]))
for file in range(len(files_albr)):
    alberio = fits.getdata(files_albr[file])
    data_alb[file] = (alberio - master_bias - (exptime_alberio/exptime_dark_40)*master_dark_40)/(master_norm_flat)

/tmp/ipykernel_5556/3063012044.py:4: RuntimeWarning: divide by zero encountered in true_divide
  data_alb[file] = (alberio - master_bias - (exptime_alberio/exptime_dark_40)*master_dark_40)/(master_norm_flat)
/tmp/ipykernel_5556/3063012044.py:4: RuntimeWarning: invalid value encountered in true_divide
  data_alb[file] = (alberio - master_bias - (exptime_alberio/exptime_dark_40)*master_dark_40)/(master_norm_flat)
In [82]:
final_alb = np.median(data_alb, axis=0)*gain/exptime_alberio
imshow_norm(final_alb[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')

/home/savannahgramze/anaconda3/lib/python3.9/site-packages/numpy/core/_methods.py:179: RuntimeWarning: invalid value encountered in reduce
  ret = umr_sum(arr, axis, dtype, out, keepdims, where=where)
Out[82]:
(<matplotlib.image.AxesImage at 0x7f2aa4093430>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ac2211d00>)
In [83]:
exptime_deneb = fits.getheader(files_denb[0])['EXPTIME']
files_denb
Out[83]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/deneb/deneb_40s_1.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/deneb/deneb_40s_2.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/deneb/deneb_40s_3.FIT']
In [84]:
data_denb = np.zeros((len(files_denb), ccd_shape[0], ccd_shape[1]))
for file in range(len(files_denb)):
    deneb = fits.getdata(files_denb[file])
    data_denb[file] = (deneb - master_bias - (exptime_deneb/exptime_dark_40)*master_dark_40)/(master_norm_flat)

/tmp/ipykernel_5556/631595824.py:4: RuntimeWarning: divide by zero encountered in true_divide
  data_denb[file] = (deneb - master_bias - (exptime_deneb/exptime_dark_40)*master_dark_40)/(master_norm_flat)
/tmp/ipykernel_5556/631595824.py:4: RuntimeWarning: invalid value encountered in true_divide
  data_denb[file] = (deneb - master_bias - (exptime_deneb/exptime_dark_40)*master_dark_40)/(master_norm_flat)
In [85]:
final_denb = np.median(data_denb, axis=0)*gain/exptime_deneb
imshow_norm(final_denb[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[85]:
(<matplotlib.image.AxesImage at 0x7f2aa4008e80>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2aa408abb0>)
In [86]:
exptime_epsP = fits.getheader(files_epsP[0])['EXPTIME']
files_epsP
Out[86]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/epsPeg/eps_Peg_40s.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/epsPeg/eps_Peg_40s_2.FIT']
In [87]:
data_epsP = np.zeros((len(files_epsP), ccd_shape[0], ccd_shape[1]))
for file in range(len(files_epsP)):
    epsP = fits.getdata(files_epsP[file])
    data_epsP[file] = (epsP - master_bias - (exptime_epsP/exptime_dark_40)*master_dark_40)/(master_norm_flat)

/tmp/ipykernel_5556/814020357.py:4: RuntimeWarning: divide by zero encountered in true_divide
  data_epsP[file] = (epsP - master_bias - (exptime_epsP/exptime_dark_40)*master_dark_40)/(master_norm_flat)
/tmp/ipykernel_5556/814020357.py:4: RuntimeWarning: invalid value encountered in true_divide
  data_epsP[file] = (epsP - master_bias - (exptime_epsP/exptime_dark_40)*master_dark_40)/(master_norm_flat)
In [88]:
final_epsP = np.median(data_epsP, axis=0)*gain/exptime_epsP
imshow_norm(final_epsP[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[88]:
(<matplotlib.image.AxesImage at 0x7f2a6491c910>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2aa40f82e0>)
In [89]:
exptime_euro = fits.getheader(files_euro[0])['EXPTIME']
files_euro
Out[89]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/europa/europa_40s_1.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/europa/europa_40s_2.FIT']
In [90]:
data_euro = np.zeros((len(files_euro), ccd_shape[0], ccd_shape[1]))
for file in range(len(files_euro)):
    euro = fits.getdata(files_euro[file])
    data_euro[file] = (euro - master_bias - (exptime_euro/exptime_dark_40)*master_dark_40)/(master_norm_flat)

/tmp/ipykernel_5556/2839034961.py:4: RuntimeWarning: divide by zero encountered in true_divide
  data_euro[file] = (euro - master_bias - (exptime_euro/exptime_dark_40)*master_dark_40)/(master_norm_flat)
/tmp/ipykernel_5556/2839034961.py:4: RuntimeWarning: invalid value encountered in true_divide
  data_euro[file] = (euro - master_bias - (exptime_euro/exptime_dark_40)*master_dark_40)/(master_norm_flat)
In [91]:
final_euro = np.median(data_euro, axis=0)*gain/exptime_euro
imshow_norm(final_euro[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[91]:
(<matplotlib.image.AxesImage at 0x7f2a6487e460>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2a64930f10>)
In [92]:
exptime_jupi = fits.getheader(files_jupi[0])['EXPTIME']
files_jupi
Out[92]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/jupiter/Jupiter_30s_1.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/jupiter/Jupiter_30s_2.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/jupiter/Jupiter_30s_3.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/jupiter/jupiter_30s_4.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/jupiter/jupiter_30s_5.FIT']
In [93]:
data_jupi = np.zeros((len(files_jupi), ccd_shape[0], ccd_shape[1]))
for file in range(len(files_jupi)):
    jupi = fits.getdata(files_jupi[file])
    data_jupi[file] = (jupi - master_bias - (exptime_jupi/exptime_dark_40)*master_dark_40)/(master_norm_flat)

/tmp/ipykernel_5556/2110104804.py:4: RuntimeWarning: divide by zero encountered in true_divide
  data_jupi[file] = (jupi - master_bias - (exptime_jupi/exptime_dark_40)*master_dark_40)/(master_norm_flat)
/tmp/ipykernel_5556/2110104804.py:4: RuntimeWarning: invalid value encountered in true_divide
  data_jupi[file] = (jupi - master_bias - (exptime_jupi/exptime_dark_40)*master_dark_40)/(master_norm_flat)
In [94]:
final_jupi = np.median(data_jupi, axis=0)*gain/exptime_jupi
imshow_norm(final_jupi[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[94]:
(<matplotlib.image.AxesImage at 0x7f2a60ccdd30>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2a648e78b0>)
In [95]:
exptime_sunn = fits.getheader(files_sunn[0])['EXPTIME']
files_sunn
Out[95]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/sun/sun_20s001.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/sun/sun_20s002.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/sun/sun_20s003.FIT']
In [96]:
data_sunn = np.zeros((len(files_sunn), ccd_shape[0], ccd_shape[1]))
for file in range(len(files_sunn)):
    sunn = fits.getdata(files_sunn[file])
    data_sunn[file] = (sunn - master_bias - (exptime_sunn/exptime_dark_40)*master_dark_40)/(master_norm_flat)

/tmp/ipykernel_5556/2247231135.py:4: RuntimeWarning: divide by zero encountered in true_divide
  data_sunn[file] = (sunn - master_bias - (exptime_sunn/exptime_dark_40)*master_dark_40)/(master_norm_flat)
/tmp/ipykernel_5556/2247231135.py:4: RuntimeWarning: invalid value encountered in true_divide
  data_sunn[file] = (sunn - master_bias - (exptime_sunn/exptime_dark_40)*master_dark_40)/(master_norm_flat)
In [97]:
final_sunn = np.median(data_sunn, axis=0)*gain/exptime_sunn
imshow_norm(final_sunn[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[97]:
(<matplotlib.image.AxesImage at 0x7f2a60cac8e0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2a60c79d90>)
In [98]:
exptime_vega = fits.getheader(files_vega[0])['EXPTIME']
files_vega
Out[98]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/vega/vega_1.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabSixData/othergroup/vega/vega_2.FIT']
In [99]:
data_vega = np.zeros((len(files_vega), ccd_shape[0], ccd_shape[1]))
for file in range(len(files_vega)):
    vega = fits.getdata(files_vega[file])
    data_vega[file] = (vega - master_bias - (exptime_vega/exptime_dark_30)*master_dark_30)/(master_norm_flat)

/tmp/ipykernel_5556/2204874859.py:4: RuntimeWarning: divide by zero encountered in true_divide
  data_vega[file] = (vega - master_bias - (exptime_vega/exptime_dark_30)*master_dark_30)/(master_norm_flat)
/tmp/ipykernel_5556/2204874859.py:4: RuntimeWarning: invalid value encountered in true_divide
  data_vega[file] = (vega - master_bias - (exptime_vega/exptime_dark_30)*master_dark_30)/(master_norm_flat)
In [100]:
final_vega = np.median(data_vega, axis=0)*gain/exptime_vega
imshow_norm(final_vega[375:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[100]:
(<matplotlib.image.AxesImage at 0x7f2a5e85b2b0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2a60c797f0>)
In [101]:
gain = fits.getheader(files_vega[0])['EGAIN']

Spectra¶

In [ ]:
In [102]:
em_alb = np.nanmedian(final_alb[370:425,:], axis=0)
fig = plt.figure(figsize=(12,2))
imshow_norm(final_alb[370:430,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[102]:
(<matplotlib.image.AxesImage at 0x7f2a5e832df0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2a5e7ffc70>)
In [103]:
fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em_alb, color='black')
plt.axvline(px_H_633, color='tab:blue')
plt.axvline(px_H_382, color='tab:blue')
#plt.axvline(H_balmer62, color='tab:orange')
#plt.axvline(H_balmer72, color='tab:orange')
#plt.axvline(H_balmer82, color='tab:orange')
plt.axvline(He_line, color='tab:red')

plt.title('Spectrum of Alberieo A')
plt.ylabel('Flux (e-/s)')
plt.xlabel('Wavelength (nm)')
plt.text(px_H_633+3, -.4, 'Balmer (5-2)', rotation='vertical', color='tab:blue')
plt.text(px_H_382+3, -.4, 'Balmer (4-2)', rotation='vertical', color='tab:blue')
#plt.text(H_balmer62+3, -2, 'Balmer (6-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer72+3, -2, 'Balmer (7-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer82-10, -2, 'Balmer (8-2)', rotation='vertical', color='tab:orange')
plt.text(He_line+3, -.4, 'He I', rotation='vertical', color='tab:red')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/alberioA_spectrum.pdf', bbox_inches='tight')
In [104]:
em_denb = np.nanmedian(final_denb[370:425,:], axis=0)
fig = plt.figure(figsize=(12,2))
imshow_norm(final_denb[370:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[104]:
(<matplotlib.image.AxesImage at 0x7f2aa41ef0a0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2a5e7dfca0>)
In [105]:
fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em_denb, color='black')
plt.axvline(px_H_633, color='tab:blue')
plt.axvline(px_H_382, color='tab:blue')
#plt.axvline(H_balmer62, color='tab:orange')
#plt.axvline(H_balmer72, color='tab:orange')
#plt.axvline(H_balmer82, color='tab:orange')
plt.axvline(He_line, color='tab:red')

plt.title('Spectrum of Deneb')
plt.ylabel('Flux (e-/s)')
plt.xlabel('Wavelength (nm)')
plt.text(px_H_633+3, -.4, 'Balmer (5-2)', rotation='vertical', color='tab:blue')
plt.text(px_H_382+3, -.4, 'Balmer (4-2)', rotation='vertical', color='tab:blue')
#plt.text(H_balmer62+3, -2, 'Balmer (6-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer72+3, -2, 'Balmer (7-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer82-10, -2, 'Balmer (8-2)', rotation='vertical', color='tab:orange')
plt.text(He_line+3, -.4, 'He I', rotation='vertical', color='tab:red')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/deneb_spectrum.pdf', bbox_inches='tight')
In [106]:
em_epsP = np.nanmedian(final_epsP[370:425,:], axis=0)
fig = plt.figure(figsize=(12,2))
imshow_norm(final_epsP[370:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[106]:
(<matplotlib.image.AxesImage at 0x7f2aafc09430>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2aa410a460>)
In [107]:
fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em_epsP, color='black')
plt.axvline(px_H_633, color='tab:blue')
plt.axvline(px_H_382, color='tab:blue')
#plt.axvline(H_balmer62, color='tab:orange')
#plt.axvline(H_balmer72, color='tab:orange')
#plt.axvline(H_balmer82, color='tab:orange')
plt.axvline(He_line, color='tab:red')

plt.title('Spectrum of Epsilon Peg')
plt.ylabel('Flux (e-/s)')
plt.xlabel('Wavelength (nm)')
plt.text(px_H_633+3, -1, 'Balmer (5-2)', rotation='vertical', color='tab:blue')
plt.text(px_H_382+3, -1, 'Balmer (4-2)', rotation='vertical', color='tab:blue')
#plt.text(H_balmer62+3, -2, 'Balmer (6-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer72+3, -2, 'Balmer (7-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer82-10, -2, 'Balmer (8-2)', rotation='vertical', color='tab:orange')
plt.text(He_line+3, -1, 'He I', rotation='vertical', color='tab:red')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/epsPeg_spectrum.pdf', bbox_inches='tight')
In [108]:
em_euro = np.nanmedian(final_euro[370:425,:], axis=0)
fig = plt.figure(figsize=(12,2))
imshow_norm(final_euro[370:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[108]:
(<matplotlib.image.AxesImage at 0x7f2aa4198d30>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ac222ce20>)
In [109]:
fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em_euro, color='black')
plt.axvline(px_H_633, color='tab:blue')
plt.axvline(px_H_382, color='tab:blue')
#plt.axvline(H_balmer62, color='tab:orange')
#plt.axvline(H_balmer72, color='tab:orange')
#plt.axvline(H_balmer82, color='tab:orange')
#plt.axvline(He_line, color='tab:red')

plt.title('Spectrum of Europa')
plt.ylabel('Flux (e-/s)')
plt.xlabel('Wavelength (nm)')
plt.text(px_H_633+3, -.3, 'Balmer (5-2)', rotation='vertical', color='tab:blue')
plt.text(px_H_382+3, -.3, 'Balmer (4-2)', rotation='vertical', color='tab:blue')
#plt.text(H_balmer62+3, -2, 'Balmer (6-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer72+3, -2, 'Balmer (7-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer82-10, -2, 'Balmer (8-2)', rotation='vertical', color='tab:orange')
#plt.text(He_line+3, -2, 'He I', rotation='vertical', color='tab:red')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/europa_spectrum.pdf', bbox_inches='tight')
In [110]:
em_jupi = np.nanmedian(final_jupi[370:425,:], axis=0)
fig = plt.figure(figsize=(12,2))
imshow_norm(final_jupi[370:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[110]:
(<matplotlib.image.AxesImage at 0x7f2a5e6886a0>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2a5e7ab1c0>)
In [111]:
fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em_jupi, color='black')
plt.axvline(px_H_633, color='tab:blue')
plt.axvline(px_H_382, color='tab:blue')
#plt.axvline(H_balmer62, color='tab:orange')
#plt.axvline(H_balmer72, color='tab:orange')
#plt.axvline(H_balmer82, color='tab:orange')
plt.axvline(He_line, color='tab:red')

plt.title('Spectrum of Jupiter')
plt.ylabel('Flux (e-/s)')
plt.xlabel('Wavelength (nm)')
plt.text(px_H_633+3, -2, 'Balmer (5-2)', rotation='vertical', color='tab:blue')
plt.text(px_H_382+3, -2, 'Balmer (4-2)', rotation='vertical', color='tab:blue')
#plt.text(H_balmer62+3, -2, 'Balmer (6-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer72+3, -2, 'Balmer (7-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer82-10, -2, 'Balmer (8-2)', rotation='vertical', color='tab:orange')
plt.text(He_line+3, -2, 'He I', rotation='vertical', color='tab:red')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/jupiter_spectrum.pdf', bbox_inches='tight')
In [112]:
em_sunn = np.nanmedian(final_sunn[370:425,:], axis=0)
fig = plt.figure(figsize=(12,2))
imshow_norm(final_sunn[370:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[112]:
(<matplotlib.image.AxesImage at 0x7f2a5e60b100>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2a5e6b9700>)
In [113]:
fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em_sunn, color='black')
plt.axvline(px_H_633, color='tab:blue')
plt.axvline(px_H_382, color='tab:blue')
plt.axvline(H_balmer62, color='tab:blue')
#plt.axvline(H_balmer72, color='tab:orange')
#plt.axvline(H_balmer82, color='tab:orange')
plt.axvline(He_line, color='tab:red')

plt.title('Spectrum of the Sun')
plt.ylabel('Flux (e-/s)')
plt.xlabel('Wavelength (nm)')
plt.text(px_H_633+3, -2, 'Balmer (5-2)', rotation='vertical', color='tab:blue')
plt.text(px_H_382+3, -2, 'Balmer (4-2)', rotation='vertical', color='tab:blue')
plt.text(H_balmer62+3, -2, 'Balmer (6-2)', rotation='vertical', color='tab:blue')
#plt.text(H_balmer72+3, -2, 'Balmer (7-2)', rotation='vertical', color='tab:orange')
#plt.text(H_balmer82-10, -2, 'Balmer (8-2)', rotation='vertical', color='tab:orange')
plt.text(He_line+3, -2, 'He I', rotation='vertical', color='tab:red')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/solar_spectrum.pdf', bbox_inches='tight')
In [114]:
em_vega = np.nanmedian(final_vega[370:425,:], axis=0)
fig = plt.figure(figsize=(12,2))
imshow_norm(final_vega[370:425,:], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/raw_vega_spectrum.pdf', bbox_inches='tight')
In [115]:
fig = plt.figure(figsize=(9,4))
plt.plot(wavelengths, em_vega, color='black')
plt.axvline(px_H_633, color='tab:blue')
plt.axvline(px_H_382, color='tab:blue')
plt.axvline(H_balmer62, color='tab:blue')
plt.axvline(H_balmer72, color='tab:blue')
plt.axvline(H_balmer82, color='tab:blue')
plt.axvline(He_line, color='tab:red')

plt.title('Spectrum of Vega')
plt.ylabel('Flux (e-/s)')
plt.xlabel('Wavelength (nm)')
plt.text(px_H_633+3, -2, 'Balmer (5-2)', rotation='vertical', color='tab:blue')
plt.text(px_H_382+3, -2, 'Balmer (4-2)', rotation='vertical', color='tab:blue')
plt.text(H_balmer62+3, -2, 'Balmer (6-2)', rotation='vertical', color='tab:blue')
plt.text(H_balmer72+3, -2, 'Balmer (7-2)', rotation='vertical', color='tab:blue')
plt.text(H_balmer82-8, -2, 'Balmer (8-2)', rotation='vertical', color='tab:blue')
plt.text(He_line+3, -2, 'He I', rotation='vertical', color='tab:red')
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/vega_spectrum.pdf', bbox_inches='tight')
In [116]:
pwd
Out[116]:
'/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix'
In [117]:
#fig, axs = plt.subplots(7, sharex=True)

fig = plt.figure()
gs = fig.add_gridspec(7, hspace=0)
axs = gs.subplots(sharex=True, sharey=False)
fig.set_figheight(8)
fig.set_figwidth(10)
axs[0].set_xlim(left=np.min(wavelengths), right=np.max(wavelengths))

axs[0].plot(wavelengths, em_sunn, color='black', label='Sun G2V')
axs[1].plot(wavelengths, em_vega, color='black', label='Vega A0V')
axs[2].plot(wavelengths, em_denb, color='black', label='Deneb A2Ia')
axs[3].plot(wavelengths, em_alb,  color='black', label='Albireo K2II')
axs[4].plot(wavelengths, em_epsP, color='black', label='Epsilon Peg K2 Ib')
axs[5].plot(wavelengths, em_jupi, color='black', label='Jupiter Gas Giant Planet')
axs[6].plot(wavelengths, em_euro, color='black', label='Europa Terrestrial Moon')

axs[1].text(px_H_633+3,   -1.5, 'H $\gamma$ (5-2)', rotation='vertical', color='tab:blue')
axs[1].text(px_H_382+3,   -1.5, 'H β (4-2)', rotation='vertical', color='tab:blue')
axs[1].text(H_balmer62+3, -1.5, 'H $\delta$ (6-2)', rotation='vertical', color='tab:blue')
axs[1].text(H_balmer72+2, -1.5, 'H ε (7-2)', rotation='vertical', color='tab:blue')
axs[1].text(H_balmer82-6, -1.5, 'H ζ (8-2)', rotation='vertical', color='tab:blue')
axs[1].text(He_line+3,    -1.5, 'He I', rotation='vertical', color='tab:red')

for ax in fig.get_axes():
    ax.label_outer()
    ax.legend()
    ax.axvline(px_H_633, color='tab:blue')
    ax.axvline(px_H_382, color='tab:blue')
    ax.axvline(H_balmer62, color='tab:blue')
    ax.axvline(H_balmer72, color='tab:blue')
    ax.axvline(H_balmer82, color='tab:blue')
    ax.axvline(He_line, color='tab:red')
    
fig.supylabel('Flux (e-/s/px)')
fig.supxlabel('Wavelength (nm)')    
    
    
fig.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/plots/science_spectrum.pdf', bbox_inches='tight')
In [118]:
#ImageNormalize(final_sunn[370:425,:], interval=ZScaleInterval(), stretch=viz.LinearStretch())
In [119]:
fig = plt.figure()
gs = fig.add_gridspec(7, hspace=0, wspace=0)
axs =gs.subplots(sharex=True, sharey=True)
fig.set_figheight(8)
fig.set_figwidth(10)
#axs[0].set_xlim(left=np.min(wavelengths), right=np.max(wavelengths))

#axs[0].imshow(final_alb[375:425])
axs[0].imshow(final_sunn[370:425,:], norm=ImageNormalize(final_sunn[370:425,:], interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
axs[1].imshow(final_vega[370:425,:], norm=ImageNormalize(final_vega[370:425,:], interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
axs[2].imshow(final_denb[370:425,:], norm=ImageNormalize(final_denb[370:425,:], interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
axs[3].imshow(final_alb[370:425,:] , norm=ImageNormalize(final_alb[370:425,:] , interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
axs[4].imshow(final_epsP[370:425,:], norm=ImageNormalize(final_epsP[370:425,:], interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
axs[5].imshow(final_jupi[370:425,:], norm=ImageNormalize(final_jupi[370:425,:], interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
axs[6].imshow(final_euro[370:425,:], norm=ImageNormalize(final_euro[370:425,:], interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')

for ax in fig.get_axes():
    ax.label_outer()
fig.tight_layout()
In [120]:
temp_sunn = 5772*u.K
temp_vega = 9602*u.K
temp_denb = 8525*u.K
temp_alb  = 4383*u.K
temp_epsP = 3963*u.K
temp_jupi = 165*u.K
temp_euro = 110*u.K
In [121]:
wein = (2.898*10**(-3)*u.m*u.K).to(u.nm*u.K)
In [122]:
peak_sunn = wein/temp_sunn
peak_vega = wein/temp_vega
peak_denb = wein/temp_denb
peak_alb  = wein/temp_alb 
peak_epsP = wein/temp_epsP
peak_jupi = wein/temp_jupi
peak_euro = wein/temp_euro
In [123]:
peak_sunn
Out[123]:
$502.079 \; \mathrm{nm}$
In [124]:
peak_vega
Out[124]:
$301.81212 \; \mathrm{nm}$
In [125]:
peak_denb
Out[125]:
$339.94135 \; \mathrm{nm}$
In [126]:
peak_alb
Out[126]:
$661.19097 \; \mathrm{nm}$
In [127]:
peak_epsP
Out[127]:
$731.26419 \; \mathrm{nm}$
In [128]:
peak_jupi
Out[128]:
$17563.636 \; \mathrm{nm}$
In [129]:
peak_euro
Out[129]:
$26345.455 \; \mathrm{nm}$

Observation Log¶

In [130]:
files = [files_bias[0],
            files_dark_07[0],
            files_dark_15[0],
            files_dark_30[0],
            files_dark_40[0],
            files_flat[0],
            files_cali[0],
            files_cali[1],
            files_albr[0],
            files_denb[0],
            files_epsP[0],
            files_euro[0],
            files_jupi[0],
            files_sunn[0],
            files_vega[0],
            files_bias2[0],
            files_dark2_10[0],
            files_dark2_20[0],
            files_dark2_40[0],
            files_dark2_50[0],
            files_dark2_60[0],   
            files_flat2[0],
            files_cali2[0],
            files_cali2[1],
            files_cali2[2],
        ]
target = ['Bias', 'Dark', 'Dark', 'Dark', 'Dark', 'Flat', 'Hg', 'Ne', 'Albireo', 'Deneb', 'Epsilon Peg', 
          'Europa', 'Jupiter', 'Moon/Sun', 'Vega', 'Bias', 'Dark', 'Dark', 'Dark', 'Dark', 'Dark', 'Flat', 'H', 'Hg', 'Ne']
exp_time = []
time_obs = []
num_expo = [len(files_bias),
            len(files_dark_07),
            len(files_dark_15),
            len(files_dark_30),
            len(files_dark_40),
            len(files_flat),
            1,#files_cali[0],
            1,#files_cali[1],
            len(files_albr),
            len(files_denb),
            len(files_epsP),
            len(files_euro),
            len(files_jupi),
            len(files_sunn),
            len(files_vega),
            len(files_bias2),
            len(files_dark2_10),
            len(files_dark2_20),
            len(files_dark2_40),
            len(files_dark2_50),
            len(files_dark2_60),  
            len(files_flat2),
            1,#files_cali2[0],
            1,#files_cali2[1],
            1#files_cali2[2],
           ]
ra =  ['-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-','-','-']
dec = ['-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-', '-','-','-']
names = ['Target', 'Exposure Time', 'Obs. Time', 'Number of Exposures', 'Right Ascension', 'Declination']

for file in files:
    head = fits.getheader(file)
    exp_time.append(head['EXPTIME'])
    time_obs.append(head['DATE-OBS'])
    #num_expo.append(len(file))

#tblr_obs = Table()
In [131]:
# Albireo
ra[8]  = '19h 30m 43.286s'
dec[8] = '+27$^{\circ}$ 57\' 34.84\"'

# Deneb
ra[9]  = '20h 41m 25.9s'
dec[9] = '+45$^{\circ}$ 16\' 49\"'

# Eps Peg
ra[10]  = '21h 44m 11.15614s'
dec[10] = '+09$^{\circ}$ 52\' 30.0311\"'

# Vega
ra[14]  = '18h 36m 56.33635s'
dec[14] = '+38$^{\circ}$ 47\' 01.2802\"'
In [132]:
jupi_timeobs = Time(fits.getheader(files_jupi[0])['DATE-OBS'])
moon_timeobs = Time(fits.getheader(files_sunn[0])['DATE-OBS'])
lat = 29.651634*u.deg
lon = -82.324829*u.deg
loc = EarthLocation(lat=lat, lon=lon, height=117*u.m)
In [133]:
jupi_body = get_body('jupiter', jupi_timeobs, loc)
In [134]:
jupi_body.fk5.ra
Out[134]:
$8^\circ57{}^\prime27.45889915{}^{\prime\prime}$
In [135]:
jupi_body.fk5.dec
Out[135]:
$2^\circ26{}^\prime44.22125532{}^{\prime\prime}$
In [136]:
ra[12]  = '8$^{\circ}$57\'27.45889915\"'#jupi_body.fk5.ra
dec[12] = '2$^{\circ}$26\'44.22125532\"'#jupi_body.fk5.dec
In [137]:
moon_body = get_body('moon', moon_timeobs, loc)
In [138]:
moon_body.fk5.ra
Out[138]:
$64^\circ55{}^\prime00.7924658{}^{\prime\prime}$
In [139]:
moon_body.fk5.dec
Out[139]:
$21^\circ26{}^\prime04.4588829{}^{\prime\prime}$
In [140]:
ra[13]  = '64$^{\circ}$55\'00.7924658\"' #moon_body.fk4.ra
dec[13] = '21$^{\circ}$26\'04.4588829\"' #moon_body.fk4.dec
In [141]:
tblr_obs = Table([target, exp_time*u.second, time_obs, num_expo, ra, dec], names=names)
In [144]:
tblr_obs.write('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabSix/tblr_obs.txt', format='latex', overwrite=True)
In [ ]: